Equilibrium glassy phase in a polydisperse hard sphere system 
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The phase diagram of a polydisperse hard sphere system is examined by numerical minimization 
of a discretized form of the Ramakrishnan-Yussouff free energy functional. Crystalline and glassy 
local minima of the free energy are located and the phase diagram in the density-polydispersity 
plane is mapped out by comparing the free energies of different local minima. The crystalline 
phase disappears and the glass becomes the equilibrium phase beyond a "terminal" value of the 
polydispersity. A crystal to glass transition is also observed as the density is increased at high 
polydispersity. The phase diagram obtained in our study is qualitatively similar to that of hard 
spheres in a quenched random potential. 
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Colloidal suspensions of polystyrene spheres coated 
with thin polymeric layers effectively behave as hard 
sphere systems jlj, exhibiting fluid and crystalline phases 
in equilibrium. For such systems, glass (amorphous 
solid) is believed to be a metastable phase 0, result- 
ing from structural arrest. However, in recent theoret- 
ical studies 0, Q of hard spheres in the presence of a 
quenched random potential, an equilibrium glassy phase 
was observed at high disorder strengths. While this phe- 
nomenon has not yet been confirmed experimentally, it 
brings forth the question of whether even the presence of 
annealed disorder in a system of hard spheres can result 
in an equilibrium glassy phase. The annealed disorder 
can be realized if the hard spheres have different sizes, 
which is the case for most colloidal suspensions. The size 
dispersity can be modeled by assuming that the diame- 
ters (cr) are sampled from a continuous distribution p(cr) , 
characterized by a parameter S, known as the polydisper- 
sity, and defined as the ratio of the standard deviation 
and the mean of the diameter distribution. In this pa- 
per we present calculations that show that for a system 
of polydisperse hard spheres, one does obtain an equi- 
librium glassy phase. In fact, our results suggest that 
for hard spheres, quenched and annealed disorder lead to 
qualitatively similar phase behavior. 

The equilibrium phase behavior of polydisperse hard 
spheres have generated a lot of interest in recent times. 
Experiments have shown that for such a system, crystal- 
lization does not occur beyond a terminal polydispersity 
of St « 0.120. Particle simulations and density func- 
tional studies|3 indicate a similar behavior, except that 
there is no consensus on the value of St ■ 

It has been observed numerically that the typical 
height of nucleation barriers for the formation of poly- 
disperse crystals grows anomalously as the polydisper- 
sity increases, thus suppressing crystallization at high 
polydispersities H . This feature has been utilized in 
experiments j2l LHj by using polydisperse hard spheres 
as a suitable system for studying the glass transition. 
Molecular dynamics simulations [llj of a supercooled sys- 



tem of polydisperse hard spheres also suggest the pres- 
ence of dynamical heterogeneity, believed to be a signa- 
ture of glassiness. 

Recently, free energy calculations found the occur- 
rence of re-entrant melting (transition from crystal to 
liquid as the density is increased) near S = St- This 
conclusion is opposed by calculations that suggest that 
the equilibrium phase at high polydispersity corresponds 
to a fractionated state|l3j. Till now, numerical evidence 
for fractionation has been found only in simulations^) 
in the grand canonical ensemble, which is different from 
the situation in typical experiments. In a recent experi- 
ment on poly disperse hard-sphere colloids in a confined 
geometry re-entrant melting was observed at the col- 
loidal monolayer adjacent to the surface, although no 
crystallization occurred in the bulk. It was claimed that 
such a phenomena could be observed in the bulk as well, 
only if one waited for sufficiently long time to allow crys- 
tallization to occur. The presence of the wall resulted 
in lowering of the nucleation barrier, thereby allowing 
crystallization and re-entrant melting to occur. 

Although there have been suggestions 0, 0] that at 
high densities or at high polydispersities, the equilibrium 
phase is a glass, there have been no calculations to sub- 
stantiate this possibility. We have addressed this ques- 
tion within the framework of density functional theory 
(DFT), using the Ramakrishnan-Yussouff (RY) free en- 
ergy functional jig. In DFT, the free energy is expressed 
as a functional of the time-averaged local density p(r). 
The RY free energy functional is given by : 

PF = J dv{p(v)\n{p{v)/ Po )-5p(v)} 

\JdvJ dv'C{\v-v'\)Spiv)Sp{v'). (1) 

Here, we have defined Sp(r) = p(r) — p as the deviation 
of p(r) from po, the density of the uniform liquid, and 
taken the zero of the free energy at its uniform liquid 
value. In EqQJ, — l/(ksT), T is the temperature and 
C(r) is the direct pair correlation function of the uniform 
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liquid at density po. 

In the polydisperse limit, the RY functional becomes^] 



pF = J dr{p(v)Hp(r)/ Po )~Sp(v)} 

- \j drj dv'C{\v-v'\)8p{v)8p{v') (2) 



where C(r) = J dat J dajp(ai)p(<jj)Cij{r), Cy-(|r - r'|) 
being the partial direct correlation function between 
spheres of size cr^ and <jj . 

In order to carry out numerical work, we discretize our 
system. We introduce for this purpose a simple cubic 
computational mesh of size L 3 with periodic boundary 
conditions. On the sites of this mesh, we define density 
variables pi = p(ri)h 3 , where p(r,) is the density at site 
i and h the spacing of the computational mesh. In terms 
of these quantities, the discretized form for the RY func- 
tional takes the form 

i 

- \'Y^2 C ij(Pi - Pl)(p 3 ~Pl), (3) 

i 3 

where the sums are over all the sites of the computa- 
tional mesh, pl = poh 3 , and Cy is the discretized form 
of the direct pair correlation function C(r) of the uni- 
form liquid. Our algor 

ithm[3 

for the minimization of 
the discretized RY functional converges to the local free- 
energy minimum whose basin of attraction contains the 
initial state. It is known that, for the monodisperse hard 
sphere system, numerical minimization of the discretized 
form of the RY functional yields a crystallization tran- 
sition at dimensionless density pi = ppa 3 ~ 0.945 if the 
mesh size h is sufficiently small |l7l Il8j. Dasgupta and 
Valls0 also obtained glassy minima for monodisperse 
hard sphere systems using this minimization scheme. 

We have extended this method of calculation to a sys- 
tem of polydisperse hard spheres. In our calculation, the 
particle size distribution is chosen to be of the Schultz 
type : p(cr) = / -f a cr a ~ 1 exp(—ja)/r(a), where a = l/S 2 
and 7 = a/ a, 8 being the polydispersity and a, the mean 
diameter of the system of particles. For the average direct 
correlation function C(r) for a polydisperse hard sphere 
system, we use the analytical expression derived using 
the Percus-Yevick approximation 20J . 

To study the stability of the crystalline minimum, we 
use the fee structure as an input for the free-energy min- 
imization. The fluid-to-crystal transition occurs when 
the free energy of the minimum obtained this way be- 
comes negative (i.e. lower than that of the uniform 
fluid). For each value of 8, the dimensionless density pi 
at which this happens is identified. The resultant phase 
diagram is shown in Figure ^ ln another method of cal- 
culating the free energyQ of the crystal, the local den- 
sity p(r) is approximated as a sum of Gaussian profiles: 



p(r) = A/7r 3 / 2 e 3 J]flexp[-(r - R) 2 /e 2 ]. In this expres- 
sion, A = (1 + ^pqVq, where rj is the density change at 
freezing and vq is the unit cell volume of a fee lattice 
with spacing a, e is the width of the Gaussians and {R} 
are the lattice points. The RY free energy functional is 
then minimized with respect to the parameters e, 77, and 
a to find the values of {/?/, 8} where the free energy of the 
crystal becomes negative. In this case, pi = p < <r 3 > is 
the polydisperse liquid density, < a 3 > being the third 
moment of the distribution p(a). As shown in Fig^ 
the results from the grid-minimization and the Gaus- 
sian approximation agree quite well, thereby establish- 
ing the correctness of the results obtained by the grid 
method. Our results corroborate earlier density func- 
tional calculations @ - there is no crystallization beyond 
a terminal polydispersity St — 0.048. Moreover, our 
calculations, using both grid-minimization and Gaus- 
sian approximation, clearly indicate that for values of 8 
slightly lower than 8t, there is re-entrant melting at high 
densities, confirming the result of an earlier free-energy 
calculation 0. 
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FIG. 1: Fluid-to-crystal transition points in the dimensionless 
density (pj)-polydispersity (8) plane. Data for two values of 
the mesh-spacing h used in the grid minimization, and the 
results obtained from the Gaussian approximation are shown. 
In all the three cases, re-entrant melting is observed at high 
density and high polydispersity. 

In the density functional calculations, the value of 
the terminal polydispersity comes out to be much lower 
than what is seen in experiments. A possible reason for 
this difference is that the average direct correlation func- 
tion C(r) used in our calculations underestimates inter- 
particle correlations for relatively large values of 8. The 
use of a more accurate C(r) would probably shift the ter- 
minal polydispersity to a higher value without changing 
other features of the phase diagram. 

The question that now needs to be answered is: What 
is the equilibrium phase of the polydisperse hard sphere 
system at high densities, beyond the point of re-entrant 
melting, and also at polydispersities where crystallization 
does not occur? Does it remain liquid or does it become 
a glass? To answer this question, one needs to locate 
glassy minima of the free energy and compare their free 
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energy with that of the uniform liquid. 

In the density functional framework, the glassy state 
corresponds to local minima with inhomogeneous de nsity 
distributions that exhibit strong non-periodicity 
To locate such minima, one first needs to construct ap- 
propriate density configurations {p^ to be used as in- 
puts for the minimization. In a recent density-functional 
study of the glass transition in monodisperse hard sphere 
systems |2l|. particle configurations from molecular dy- 
namics (MD) simulation were used to produce the den- 
sity field for calculating the free energy using the Gaus- 
sian approximation for the local density. Similarly, in 
our discretized method for doing the DFT calculation, 
the input density field was constructed from MD simula- 
tion of polydisperse hard spheres. The simulations were 
carried out for 470 particles with their diameters sampled 
from the Schultz distribution with 8 — 0.0289. Using a 
modified form of the Stillinger-Lubachevsky compression 
algorithm [2^, highly dense configurations of polydis- 
perse hard spheres could be created. The input density 
field {pi} was calculated from these configurations using 
a computational grid of mesh-spacing h w 0.055\ The 
RY functional was then minimized in the space of the re- 
sulting 4.096 x 10 6 density variables to obtain the density 
configuration at the local minimum of the free energy. As 
in the case of monodisperse hard spheres j2^|, we obtain, 
in this case also, free-energy minima with the {pi} having 
a glassy structure. The structure of a local minimum may 
be characterized by the two-point correlation function 
g(r) of the local density variables {pt] at the minimum. 
This function is defined as g(r) —< p(0)p(r)/p 2 >. In 
Fig. |2 we have plotted the g(r) for a glassy free energy 
minimum obtained for 8 — 0.0289. The glassy nature of 
the density distribution is indicated by the split second 
peak of the g(r), with the sub-peaks occurring at 1.77a 
and 2.02cr. These numbers are close to those for a dense 
random packing of monodisperse hard spheres, where the 
sub-peaks occur at around 1.7a and 2(r|24j. 



The density configuration {p^ at the local minimum 
for 8 = 0.0289 was thereafter used as input to search for 
similar glassy minima at other values of the polydisper- 
sity. The density at which the free energy of the glassy 
structure becomes lower than that of the uniform liquid 
defines the point of the liquid-to-glass transition. Fig. [3| 
shows a plot of the liquid-to-glass transition density for 
various values of 8. As shown in the plot, one obtains 
glassy minima for polydispersities higher than <5 t , the ter- 
minal value beyond which crystallization does not occur. 
Also the glass transition point shifts to higher densities 
as a function of increasing polydispersity and at large 
enough polydispersity, near 5 = 0.10, the glass transi- 
tion point approaches the random close packing limit for 
polydisperse hard spheres [EHH When the density con- 
figuration for the minimum for 8 — 0.02 is used as input 
for minimizing the RY functional for 8 = 0.0, i.e. for a 
system of monodisperse hard spheres, we obtain a free 
energy minimum whose {pi} correspond to a polycrys- 
talline structure. 




FIG. 3: Glass transition points, plotted in the dimensionless 
density (pi)-polydispersity (5) plane. At these points, the RY 
free energy of the glassy minimum becomes negative. The line 
drawn through the data points is a guide to the eye. 




FIG. 2: The two point density correlation function g(r) 
for a glassy minimum obtained for a polydispersity value of 
8 — 0.0289 and dimensionless liquid density pi — 1.11. Also 
plotted in the figure are the peaks of the pair distribution 
function (not drawn in the same scale as g(r) for the glass) 
for a fee lattice at this density. 



Having determined the densities (pi) beyond which the 
crystalline and glassy phases, respectively, have lower 
free energies compared to that of the homogeneous liq- 
uid phase, the next task is to determine the final phase 
diagram, i.e., to determine which phase (crystal, glass or 
liquid) is the thermodynamically stable one at different 
points in the (6, pi) plane. At any density, if more than 
one local minima of the free energy are present, the ther- 
modynamically stable phase corresponds to the one with 
the lowest free energy. Using this criterion, we have ob- 
tained the full phase diagram of the polydisperse hard 
sphere system, which is shown in Fig. H The inset m 
Fig. 0] shows, for example, how the crystal-to-glass tran- 
sition line is located - for each polydispersity, the free 
energies of the glassy and crystalline minima at different 
values of pi are compared and the value of pi at which the 
free energy of the glassy minimum becomes lower than 
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that of the crystalline minimum defines the crystal-to- 
glass transition point. 




FIG. 4: The overall phase diagram of a system of polydisperse 
hard spheres in the dimensionless density (pj)-polydispersity 
(5) plane showing the presence of the three phases: uniform 
liquid, crystal, and glass. The lines drawn to denote the phase 
boundaries are only meant to be guides to the eye. The in- 
set shows the free energies of the the glassy and crystalline 
minima for 8 = 0.043 (shown by the dotted line in the main 
figure). The crystal is the equilibrium phase till pi ~ 1.126, 
after which the glass becomes the equilibrium phase. 

It is important to note that the re-entrant melting 
shown in Fig. [2 is not present in the final phase diagram. 
Now, as one moves along a line of fixed polydispersity 
(for example, 5 — 0.043 as shown by the dotted line in 
Fig. 0J , one first encounters the fluid phase, then a crys- 
talline phase and at even higher densities, glass becomes 
the stable phase (as is clear from Fig.|3). The other sig- 
nificant feature of the phase diagram is that, at a fixed 
high dimensionless density (pi > 1.05), if we increase the 
polydispersity 5, the crystal undergoes a transition to the 
glassy state. We find that the crystal-glass transition line 
is below the fluid-crystal re-entrant line in Fig^ A pos- 
sible reason for this behavior is that as the polydispersity 
increases, the number of smaller particles in the system 
increases and their presence makes the glassy phase to 
be entropically favorable as compared to the crystalline 
phase. Very similar features, including the presence of an 
equilibrium glassy phase and crystal-to-glass transitions 
upon increasing density or disorder, have been found0,0] 
in the phase diagram of a system of hard spheres in a 
quenched random potential. 

Our calculations do not take into account the possibil- 
ity of having fractionated phases. We need to compare 
the free energy of the glassy phase with that of fraction- 
ated phases at high polydispersities to determine which 
one is the equilibrium phase. A proper procedure for do- 
ing such calculations in the framework of DFT is not yet 
available. However, one must remember that for many 
glass-forming liquids, phase separation is possible but 
may not be observed in experimental time scales. This 
may be the case for this system as well. 

In summary, we have shown, using density functional 



theory, that re-entrant melting occurs at high densities in 
a system of polydisperse hard spheres if only the crystal 
and liquid states are considered. However, the re-entrant 
liquid state is replaced by a glassy state in the full phase 
diagram because of the presence of glassy minima with 
lower free energy. At high polydispersity, when crystal- 
lization does not occur, our calculations show that the 
glassy state can be an equilibrium phase for this system. 
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